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Abstract 

Formulating a dust filled spherically symmetric metric utilizing the 3 + 1 for- 
malism for general relativity, we show that the metric coefBcients are com- 
pletely determined by the matter distribution throughout the spacetime. Fur- 
thermore, the metric describes both inhomogeneous dust regions and also 
vacuum regions in a single coordinate patch, thus alleviating the need for 
complicated matching schemes at the interfaces. In this way, the system is 
established as an initial-boundary value problem, which has many benefits for 
its numerical evolution. We show the dust part of the metric is equivalent to 
the class of Lemaitre-Tolman-Bondi (LTB) metrics under a coordinate trans- 
formation. In this coordinate system, shell crossing singularities (SCS) are 
exhibited as fluid shock waves, and we therefore discuss possibilities for the 
dynamical extension of shell crossings through the initial point of formation 
by borrowing methods from classical fluid dynamics. This paper fills a void in 
the present literature associated with these collapse models by fully developing 
the formalism in great detail. Furthermore, the applications provide examples 
of the benefits of the present model. 



1. Introduction 

It is well established that naked singularities arise from the gravitational 
collapse of inhomogeneous, spherically symmetric dust spheres in general 
relativity [231 El IH IIS]. One type of these are shell-crossing singularities 
(SCS), which are not perceived to violate cosmic censorship due to their 
weak nature (see e.g. [SIT and references therein). However, difficulties 
often arise in the interpretation of the physical nature of SCS due to the 
four-dimensional construction of spacetime solutions usually having little or 
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no resemblance to our intuitive perception of the universe. Therefore, we 
utihze the 3 + 1 formahsm for general relativity to model the spacetime as 
an initial value problem in a single coordinate patch. In this way the SCS 
are akin to shock waves in fluid mechanics, and therefore our perception and 
understanding of the solutions are greatly enhanced. 

The 3 + 1 formalism involves specifying data on an initial Cauchy hy- 
persurface, and subsequently evolving through time. The ten Einstein field 
equations (EFEs) are decomposed into four constraint equations which must 
be satisfied on every spacelike hypersurface, and six evolution equations. 
Furthermore, they can be supplemented by the integrability conditions for 
the EFEs, which are the Bianchi identities. The once contracted Bianchi 
identities in four dimensional spacetimes provide 16 conditions, whilst the 
twice contracted provide the remaining four pieces of information, which are 
the conservation of energy- momentum equations. 

We apply the 3 + 1 formalism for the EFEs and Bianchi identities to 
the case of a spherically symmetric sphere of inhomogeneous dust. This 
article acts to fill a void in the literature created by [9l [10] by deriving the 
formalism in full detail. Furthermore, while j9] and jlOj used more general 
fluids, analytic solutions for the equations derived therein are extremely 
difficult to find, and hence the equations are awkward to interpret. It is 
therefore pertinent to simplify the fluid in order to extract analytic solutions 
such that the full strength of the formalism can be displayed. The dust 
application we present in this article is solved analytically, and hence a 
greater understanding of the structure of the equations is gained. One of 
our main results from the dust section is that the metric coefficients are 
completely determined by the energy-momentum fields. Furthermore, to 
determine the matter, and hence the entire spacetime, only the specification 
of the matter distribution on the initial hypersurface is required. 

We show our solution under a coordinate transformation is equivalent 
to the class of Lemaitre-Tolman-Bondi (LTB) solutions |12 | I2 H [3]. Further- 
more, the line element is such that an exterior vacuum region, namely the 
Schwarzschild spacetime expressed in a generalized form of the Painleve- 
Gullstrand (PG) coordinates |18| [7]. is described using the same coordinate 
patch. In this way, the work is related to that of [T] and [B] who, in a sin- 
gle coordinate patch describe the Oppenheimer-Snyder (OS) collapse [T7j . 
i.e. a Friedmann-Robertson- Walker (FRW) interior joined to an exterior 
Schwarzschild spacetime. As the FRW spacetime belongs to the class of 
LTB spacetimes, our work is a generalization of [1] and [6]. 

The evolution of the system is determined by two, first order equations 
in the matter fields. The nature of these equations is such that multi- 
valued solutions in the mass can develop given smooth initial data. These 
are equivalent to the SCS discussed, and we analyse the initial conditions 
necessary for these to occur. Furthermore, we discuss a method for extending 
these SCS beyond their point of initial formation utilizing classical fluid 
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dynamics methods. In this way, we show when a globally naked SCS forms 
(i.e. beyond the apparent horizon) it will evolve to become a locally naked 
singularity, and eventually evolve into a Schwarzschild singularity at r = 0. 

As all aspects of the solutions we present can be solved explicitly in terms 
of the initial data, the authors believe this formulation is of great value for 
testing numerical evolution codes. In particular, the evolution of such a 
system contains many complications, including the formation of SCSs, as 
well as the interface between the two regions of the spacetime. Numerical 
codes are extremely well developed to handle shock waves from classical fluid 
dynamics, and the same methods can be employed in the relativistic regime. 
However, the evolution of the interface between the two regions may be of 
more primary concern. In particular, the advantage of having both regions 
in a single coordinate patch implies no extra work in the code is required 
at the moving boundary. One simply sets up the initial Cauchy data, with 
the matter terms going to zero at some finite radius, and the evolution will 
naturally take care of the free boundary. 

The structure of the paper is as follows. In section |2] we derive the 
general form of the LTB metric and analyse the solution. In section [3] we 
look at the initial and boundary conditions required for the specification 
of the spacetime, and also look at the formation of the apparent horizon. 
In |5] we reduce the system to the marginally bound case and derive the 
initial conditions that provide shock waves. Finally, in section [6] we look at 
a specific example of the evolution of the system with a shock wave. 

Geometrized units are employed throughout whereby G = c = 1. Greek 
indices run from ... 3 and Latin from 1 ... 3. The Einstein summation 
convention is used, index conventions follow [14 and is the line element 
for the two-sphere. 

2. The Metric 

Here, we pedagogically derive the spacetime metric to be used for the 
remainder of the article. This metric is beneficial as it describes both an 
interior dust region as well as an exterior vacuum region without the need 
for complicated matching schemes across the interface. In particular, both 
regions are represented not by two separate line elements patched together, 
but as a single, all encompassing, line element. In this way, the metric we 
derive is a generalization of the PG vacuum solution to include dust in any 
region of the spacetime desired. 

The derivation of the metric utilizes the 3 + 1 formalism for general 
relativity. This is a powerful method that allows the resulting spacetime to 
be represented as an initial value problem. In this way, initial Cauchy data is 
prescribed and the subsequent evolution depends continuously on the initial 
data. Although analytic solutions are provided here, the authors believe 
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the resulting class of spacetimes provides excellent test-beds for numerical 
schemes; in particular for gravitational collapse. 

2.1. ADM equations The decomposition of the EFEs into a 3+1 system 
is otherwise known as the ADM formalism [2]. Solving the equations in this 
form establishes a metric appropriate to an initial value problem. 

Consider a four-dimensional spacetime foliated with three-dimensional 
spacelike hypersurfaces. We denote the four-coordinates with := (t, x*), 
and a spherically symmetric line element can be reduced for dust, without 
loss of generality, to 

ds'^ = - (a^ - U'^p'^) dt^ + 2U'^l3dtdr + U'^dr'^ + r'^dO^. (1) 

where a = a{t,r) is the lapse function, U = h({t,r) and (3 = (3{t,r) is the 
radial component of the shift vector. 

The energy-momentum tensor for dust is given by 

T^u = pn^Uy, (2) 

where p is the energy-density, and is the normal vector field tangent to 
the fluid lines. We demand three properties for the normal vector, given by 
two equations 

rf^na = — 1 and n\^^Vi,n^T^ = 0, (3) 

where is the unique four-covariant derivative operator. The first equation 
in (|3| implies the vector is both timelike and normalized, and the second 
equation implies the normal vector is hypersurface forming. Utilizing the 
line element ([T| and both equations in (|3]), one can show the normal vector 
in component form is given by 

n^ = - (1,-/3,0,0). (4) 
a 

Using the normal vector, the EFEs, G^^ = SirT^^, can be decomposed 
into four constraint equations which must be satisfied on each spacelike 
hypersurface. These are the Hamiltonian constraint 

^R+lK^-AjA'^ = 16np, (5) 

where is the three-Riemann scalar, Aij and K are the trace-free part and 
trace of the extrinsic curvature respectively. There are three momentum 
constraints 

D^Ai, = '^DiK, (6) 
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where Di is the unique three-covariant derivative operator. The remaining 
six equations are the evolution equations, 

2CnK =K^ + ^Ai.jA'^ + ^^R - ^DWia, (7) 
n 2 ^ 2 a ^ ^ 

- -D.Dja + D^Dka, (8) 

where is the Lie derivative operator with respect to the normal vector. 

Equations ([5]j8| are the ADM equations, and these need to be supple- 
mented with the conservation of energy-momentum equations 

v„r"^ = 0. (9) 

These can be decomposed giving the continuity equation which, for dust 
reduces to 



E: = £„(ln/5), (10) 

and the Euler equations which, in dust, are 

pDi(lna) = 0. (11) 

As we are considering systems with non-vanishing energy density, equation 
(11) along with the form of the normal vector (|4| , is enough to show the 



lapse function is only a function of the time coordinate. Utilizing coordinate 
freedom, we set the lapse function to unity without loss of generality, 

a = l. (12) 

This physically implies that the dust particles are moving along timelike 
geodesies of the spacetime, which is a result of there being zero pressure. 

2.2. Governing equations We begin by putting the metric ([T]), into 
the ADM system of equations. Noting that the Lie derivative acting on a 
scalar function, ^, is simply 



helps to recognize patterns in the systems of equations. In particular, we 
immediately find the radial component of equation ([6]) reduces to 

CnU = 0. (14) 
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This term, and its derivatives appear readily throughout the system of equa- 



tions. Substituting ( 14 ) into the remaining equations impHes the rest of the 
system reduces to a comphcated set of differential equations in lA^ p and (5. 

Applying various combinations of the 3 + 1 EFEs results in being able 
to write the system algebraically in lA and p, leaving only derivatives of the 
shift function (5. In particular, combining ^ and ([8| we find the derivatives 
of U cancel and a relation between p and derivatives of the shift results in 

/•r 

r'^Cn(3 = 47r p{t,a) a'^da := M{t,r), (15) 

Ja=0 

where we have defined a "mass" function according to a Newtonian defi- 
nition. We note here that this function does not physically represent the 
mass of the system as a Newtonian volume element has been used. How- 
ever, this function arises in a natural way. Furthermore, we will show that 
this is precisely the mass function commonly defined in the LTB spacetimes. 
Furthermore, in a vacuum region of the spacetime where the energy den- 
sity vanishes, the mass function simply becomes a constant. We will see 
that this vacuum spacetime is given by the Schwarzschild spacetime in a 



generalization of the Painleve-Gullstrand coordinates (see section 2.5). 

Adding equation ([T]) to six times the radial component of equation ^ 
enables derivatives of U to again cancel, resulting, after some algebra, in 

We now find that the final metric function U is related to the shift, and 
therefore, implicitly related to the energy density. Thus for dust, the spher- 
ically symmetric metric in the form given by ([T| is entirely dependant on 
the initial matter field. 



Now by substituting (16) into (14) we find a second order evolution 
equation solely on the shift 

= Cn (/?' - 2rCnP) . (17) 

Summarily, the line element takes the form 

_ _dt^ , m + drf , , 



where /3 is given by a solution of equation (17) and is further constrained 
by the requirement that the signature of the spacetime be Lorentzian, /J^ — 
2rCnP > -1. 

Due to the method of deriving the metric, we note that it is conducive 
to an initial value formulation. In particular, the evolution of the system 



is governed by the second order partial differential equation (17). To solve 
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this equation one must pose boundary and initial data, which is discussed 
in section |3l 

The physical realization this method has brought is due to the manner 
in which it was derived. Many methods for solving the EFEs require simply 
solving for the full four-dimensional spacetime with minimal consideration of 
the physics until the solution has been found. It is then difficult to decipher 
the physics due to the seemingly arbitrary nature of the coordinates. By 
specifying spacelike hypersurfaces and allowing for the subsequent evolution 
of the system, we have reduced the problem to one that relates closer to the 
physical intuition grasped from everyday life. 

2.3. Gravitoelectromagnetism Another important set of equations 
which can be solved alongside the ADM equations ([5]j8| and the conserva- 
tion equations 

sions one can go to the once contracted Bianchi identities without loss of 
generality. By introducing the Weyl conformal curvature tensor, the once 
contracted Bianchi identities can be expressed in terms of the energy mo- 
mentum tensor 

V"a^.. = 8^ (v[uT^]f, + ^5^[.V,]r„"^ . (19) 

The Weyl tensor can be decomposed into two spatial, trace free, gravito- 
electromagnetic (GEM) tensors, known as the electric and magnetic confor- 
mal curvatures, Eij and Bij respectively, where 

Eij : = Qaj/^n^n^, (20) 

Bij : = leaip^Cf'^Sjn''n\ (21) 

The electric conformal curvature is a measure of the tidal forces present 
in the spacetime, analogous to the tidal tensor in Newtonian gravitation. 
The magnetic conformal curvature has no Newtonian analogue, and is in 
some sense a measure of the intrinsic rotation associated with the space- 
time [5]. Therefore, while Bij = in spherical symmetry, we note that it 
will have non-trivial contributions for spacetimes with less symmetries. For 
example, one aim of the current program of research is to generalize this 
method to include the family of Robinson- Trautman solutions which con- 
tain gravitational radiation. These spacetimes will begin with a non-zero 
magnetic curvature. However one in these scenarios that this contribution 
will be "radiated away" in the form of gravitational radiation such that 
the steady-state solution is spherically symmetric. Furthermore, a radiat- 
ing, axially-symmetric spacetime will have two components of the magnetic 
curvature. One component will be associated with the intrinsic angular mo- 
mentum of the spacetime, and will remain once the system has radiated 
to a steady-state. The other component associated with the gravitational 
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TOpT) are the remaining Bianchi identities. In four dimen- 



radiation will be radiated away in the vein associated with the Robinson- 
Trautman solutions. 

As mentioned, for spherically symmetric dust the magnetic curvature 
tensor vanishes and the only contribution to the Weyl tensor is from the 
electric curvature. The trace-free nature of the electric curvature, as well as 
the spherical symmetry imposed, implies it has the simple form 



S,-'' = diag(-2A,A,A), 



(22) 



where A = X{t, r) is the tidal forces associated with the spacetime. Decom- 
posing the Bianchi identities, given Bij = 0, yields two non-trivial equations 



Svr 



A/0, 



(23) 
(24) 



Equation (23) can be integrated by using the form of the electric curvature 
(22), and by substituting the definition for the mass (15), to give the tidal 
forces for the spacetime 



d fM 



3 dr \ r 



(25) 



We note that in a vacuum region, the mass function is constant, and the 
tidal forces therefore become 



A,, 



(26) 



which is the familiar form known for the Schwarzschild spacetime. These 
quantities will become important for the discussion of shock formation in 
section 15.11 



The other non-trivial GEM equation (24) can be calculated using equa- 
tions (15) and (22), and also substituting (25), and it simply gives an inte- 
grability check for the evolution equation which will be given for the mass 
function (see section 2.4 in particular equation (35b)). 

2.4. Equivalence with LTB We have derived the reduced field equa- 
tions (17), for the line element (18) that represents spherically symmetric 
dust. We now show that this can be transformed into LTB coordinates 
{t,R,9,^). 

Let r be a function of both the new time and new radial coordinates, 
i.e., r = r(r, it!), and also let t = t, then let 



9r(r, R) 
dr 



-f3. 



(27) 
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In this new coordinate system, the normal vector has the form n'^ = 
(1,0,0,0), and the observer is now comoving with the fluid. In these co- 
ordinates the Lie derivative of a scalar function is simply the derivative with 



respect to the new temporal coordinate, r, and the condition given by (17) 
becomes 



d_ 



dr{T, R) 
dr 



+ 2r(T, R) 



d^r{T, R) 



(28) 



This is a third order differential equation in the function r{T,R), and can 
be integrated twice to yield 



dr{T, R) 
dr 



E{R) + 



2m{R) 
r{T,R)' 



(29) 



where E(R) and m{R), known as the energy and Misner-Sharp mass [13] re- 
spectively, are functions of integration. Furthermore, rearranging the metric 
we find 



dr{T,R) 
dR 



1 + E{R) 



-dR^ + r{T, RfdD? 



(30) 



Equations (30) and (29) are the standard form of the LTB metric [12 1 121 } [3]. 

The function E{R), arbitrary up to the constraint E{R) > —1, is a 
measure of the energy of a shell at a radius R [20] . Therefore, this function 
being equal to zero is equivalent to saying the particles at infinity have zero 
kinetic energy. 

It was found during the above transformation that this energy function 
can be expressed in terms of the shift function in our coordinates 



E{R) = P^- 2r£„/3, 



(31) 



which is simply the first integral of equation (17). 



Transferring back into the original coordinates, we note that the energy 



function is a function of t and r. Furthermore, substituting equation (31) 



back through (17), we see 



CnE = 0. 



(32) 



Now, substituting the mass defined in equation (15) back into (31) gives 

2M 



E + 



from which we can further show this implies 

CnM = 0. 



(33) 



(34) 
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Finally, the above equations can be put back through the line element, and 
we find the resulting system can be summarized by the metric and two 
constraint equations 



JE + 2M/r dt^dr] 

dt'+^ ^ ^ + r^dn\ 



l+E 



dM 
~dt 
dE 
Ik ~ 



E + 



2MdM 
r dr 



E + 



2M dE 

r dr 



(35a) 
(35b) 
(35c) 



Equations ( 35b ) and ( 35c ) provide a coupled system of differential equations 



which are solved concurrently to determine the dynamics of the system. 



We note that the positive root of equation ( 33 ) has been selected as this 



represents a collapsing model. Choosing the negative root is equivalent to re- 
versing the time coordinate, and therefore gives an expanding, cosmological 
type model. 

2.5. Vacuum Regions While we have derived this solution as a dust 
problem, we note that these coordinates also describe vacuum regions. This 
is shown by simply letting the energy density vanish at some radii on the 
initial hypersurface. This implies that the mass function becomes constant, 
and thus equation (35b) is trivially satisfied. Therefore, vacuum regions 



of the system are described by the line element (35a) and equation (35c 
with M = Ms G R. We describe this system as the "generalized Painleve- 
Gullstrand" (GPG) coordinate^as they reduce to the usual PG coordinates 
jl81 [7] for the particular case of -E = 0. Evaluating the Einstein tensor for 
the GPG metric we find it vanishes as expected. This is therefore a vacuum 
solution of the field equations, and spherical symmetry along with Birkhoff's 
theorem implies this solution must be diffeomorphic to the Schwarzschild 
metric. 



However, the vacuum solution described by (35a) and (35c) is actually a 



family of solutions parametrized by the energy function. Therefore, for all 
solutions of equation (35c), there will exist a coordinate transformation to 
the Schwarzschild metric, in coordinates (T,r,9,4>). This transformation is 
given by solutions to the following coupled system of differential equations 



2M,\ dt 
dr 



--1 + E, 



2M, 



+ E. 



(36a) 
(36b) 



'^A more detailed analysis of these vacuum coordinates will be presented in a future article. 
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To show this is a vahd coordinate transformation amounts to showing 



that there exists a solution of (36), which is done by checking the integra- 



bihty conditions. That is, by differentiating equation (36a I with respect to 



the radial coordinate, and equation (36b) with respect to the temporal co- 



ordinate, one can show the partial derivatives commute providing equation 
is satisfied. Therefore, providing the energy function is a solution of 
the coordinate transformation is valid, and applying the coordinate 



35c 



35c 



transformation yields 



r J 1 - 2Ms/r 



(37) 



which is exactly the Schwarzschild metric. 



3. Initial and Boundary Conditions 

We can now establish the initial and boundary data required to solve the 
system of equations that describe the spacetime. We are required to specify 
an energy density for the initial hypersurface. Due to the formulation of the 
solution, this energy density can take any form, including vanishing for finite 
regions. For example, if one wanted to study the gravitational collapse of 
a spherical body, then one specifies an appropriate energy density out to a 
finite radius, at which point the energy density is allowed to vanish beyond 
that radii. However, one can also study the collapse of concentric shells of 
matter; for example by specifying successive Heaviside step functions for the 
initial density. The robustness of the formalism allows for the study of any 
initial configuration of spherically symmetric dust. 

Whatever form the initial energy density takes, we can evaluate the 



initial mass function through the definition given by equation (15). For 
vacuum regions, this mass function reduces to a constant. We therefore have 
the initial conditions for the mass function throughout the initial spacelike 
hypersurface. 

The data for the energy function is a little more tricky. By expressing 



equation (35b) in terms of the energy function, and substituting equation 



(15), after some algebra we find 




EiO,r< re) = I , I - - / Pi0,a)a'da. (38) 



r 







Therefore, the energy function is found by specifying two terms, namely 

and p{0,r). (39) 



dp 
di 



t=o 
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While this works for the non- vacuum regions, one stih has the freedom to 
choose the form of the energy function for the vacuum regions. However, if 
ah regions of the spacetime are to be described by a single coordinate patch, 
then one can employ continuity of the energy function across the interface of 
the matter filled and vacuum regions to determine the boundary conditions 
for the vacuum region. 



4. Apparent horizon 



No discussion of gravitational collapse is complete without contemplating 
the apparent horizon. We shall now show that the coordinate system used 
herein allows for the description of the apparent horizon in a clear and 
concise manner. 

The apparent horizon is defined as the boundary of the closure of the 
union of all trapped regions on a Cauchy surface [22]. An alternate, but 
equivalent definition is given by the surface with a vanishing expansion fac- 
tor, Q, which is defined as the divergence of a congruence of null geodesies 

G := V«fc". (40) 



Here is a null vector which is everywhere tangential to the congruence 
of radial null geodesies. To derive this null vector, we must look at the 
equations of motion describing the null geodesies in the spacetime. These 
are given by the Lagrangian, C, and the Euler-Lagrange equations. 



1 



= 



l + E 
d dC 



- I 



2M 

r 



f + 2 



2M 



dX di 



dC 



and 



= 



r 

d dC 
dX dr 



+ Eir + 



dC 

dr 



0, (41) 
(42) 



where a dot denotes differentiation with respect to the affine parameter 
along the geodesies, A. Furthermore, we have used spherical symmetry to 
imply that the apparent horizon will only depend on the radial and temporal 
coordinates, i.e. 9 = (p = 0. 

The three equations ( 4l|42 ) can now be integrated to solve for i and r 
in terms of a single constant of integration and the metric coefficients. The 
null vector := can then be expressed 



k^" 



^/l + E A+E-^/l + E 



2M 



+ ^ ,0,0 



(43) 



where the constant of integration has been scaled to unity without loss of 
generality. Taking the divergence of this quantity, as indicated by equation 
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(40), one finds the expansion factor vanishes along the surface given by the 
implicit parametric equation 



r{t) = 2M{t,r{t)). (44) 

This remarkably simple form is true for all choices of initial conditions. Fur- 
thermore, at the interface between the matter and vacuum regions, the mass 
function simply becomes the Schwarzschild mass, and the horizon reduces 
to the familiar event horizon in the Schwarzschild spacetime. 



5. Marginally Bound Solution 

For the remainder of the article we only consider the case where the 
energy function, E{R), vanishes. Known as the marginally bound cas^ 
this model provides a simpler example than the non-zero energy case. The 
marginally bound system is described by the line element and single condi- 
tion on the mass 



ds'^ = -dt^ + {^'^dt + drj + r'^dn'^, (45a) 
dM [2MdM , , 



Equation (45b) is a first order, quasi-linear partial differential equation 
which governs the dynamics of the system. The fact that this can easily 
be written in conservation form is important to the numerical analysis of 
the system. In particular, systems that can be written in conservation form 
are preferred since the convergence (if it exists) to weak solutions of the 
equations are guaranteed by the Lax-Wendroff theorem 



The general solution to (45b) is given implicitly by 




M = J^{ + tV2M , (46) 



where .7-" is a function of integration. To find the particular solution, and 
hence determine the dynamics of the system, one simply requires the input 
of an initial mass distribution. This initial data can be expressed in terms of 
the energy density, which is translated into an initial mass via the definition 



(15) 



The nature of equation (46) is such that given smooth initial data, the 
solution can evolve to be multi-valued for the mass function, violating phys- 
ical intuition regarding the behaviour of mass. These multi- valued solutions 



"E{R) < is the bound case and E{R) > is the unbound case. 
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correspond to SCS in general relativity, and are also described by mathe- 
matical shock waves in this formalism. The solution prior to the formation 
of these shocks is described by the classical solution of the differential equa- 



tion (45b) while after the formation of the shock the solution can still be 



evolved utilizing a weak solution. 



5.1. Classical Solution By analyzing the characteristics of (45b), we 



can determine the point at which the classical solution fails, and also the 
initial conditions required for this to occur at some point within the collapse 
process. 



The characteristics of equation (45b) can be represented by introducing 
a parameter ^, such that 



t = t(^), r = r(0 and M = M(C). 



(47) 



Furthermore, by defining t(0) = without loss of generality, and defining a 
new parameter s := r(0), we see M(t(0),r(0)) = M(0,s) := Mo(s). Now, 
the characteristics are given by the parametric equations 



M = Mo(s), 



.3/2 



(48a) 
(48b) 



The solution becoming multi-valued corresponds to the intersection of the 
characteristic curves, which is interpreted as when the transformation (^, s) — > 
{t, r) becomes non-invertible. That is, the determinant of the Jacobian of 
transformation vanishes. Therefore, to find conditions for the formation of 
the shock, we must solve 



J 



dt 

dr 



at 

ds 

dr 
ds 



0. 



It is trivial to show from equations (48) that 



dr 

OS 



This condition is equivalent to 



a/3 



1 



(49) 



(50) 



(51) 



(52) 



where a prime denotes differentiation with respect to s. We are only con- 
cerned with shocks that occur in the first quadrant of the (t, r) plane. There- 
fore, we must determine the conditions on Mq and Mq which imply the 
Jacobian vanishes. 
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As s = r{t = 0), we see s G [0,oo). Furthermore, p > and Mq{s) 



is defined by letting t = in equation (15), implying Mq(s) > for all s. 



Therefore, by demanding t > in (51), we see 



Mq > 0. 



(53) 



This will only not hold true when p(0, r) = for some finite range of r, as 
this will imply Mg = 0. Therefore, the characteristics emanating from radii 
for which there is a vanishing energy density will never cross, as expected. 



Now, demanding r > in equation (52 ) implies, after some manipulation 

d /Mo' 



> 0. 



(54) 



ds \ s'^ 

Thus, shocks occur in the first quadrant of the (t, r) plane if for any finite 



range of s G [0,oo), the inequality in (54) is satisfied. 



Finally, letting t = in (25) implies the "tidal force" at t = is 



X{0,s) 



s d [Mo 
3 ds 



and considering s > 0, the inequality in (54) is equivalent to 



A(0,r) < 0. 

Summarily, we have proved the following theorem: 



(55) 



(56) 



Shell-Crossing theorem 5.1. SCS occur in marginally bound LTB 
collapse if and only if A(0, r) < for some finite range of r G [0, oo) . 

Conversely, if A(0, r) > for all r, then shock waves will not form. 
Furthermore, the initial time for the shock to begin, denoted tg-, can be 



determined solely from the initial conditions by taking the minimum of (51 ), 



V2 



mm 



(57) 



5.2. Weak Solution The solution is satisfied classically until the point 
where the characteristics first cross. This is the initial point of the shock 



surface, and beyond this point the solution to equation (45b) is only given 



by a weak solution. The analysis of the weak solution is made simpler by 
putting the equation into conservation form. This is done by rescaling the 
radial coordinate such that r := r^/^, implying 



(58) 



where f{M) := — \/2 Af The weak solution is defined as follows: consider 
a test function ip £ which vanishes everywhere outside a rectangle defined 



15 



by < t < r and a < r < b, and also on the lines t = T, a = r and b = r 
(for some a, b and T G R). A weak solution is defined to be a solution of 



T rb 
t=0 Jf=a 



b 



+ / Mo(f)^A(0,f)df = 0, (59) 

J f=a 

for all test functions ip. 

The weak solution contains a world-line whereby the solution contains 
a jump discontinuity which is the shock surface. This world-line can be 
expressed parametrically in terms of the time coordinate, i.e. rs(t). Either 
side of the shock, the solution is satisfied by the usual classical solution 
discussed in section 15.11 

A method for evolving shocks beyond the initial point of formation, 
utilized in many classical scenarios, is to introduce a viscosity term into the 
differential equation, 

dM d ,,,,, d'^M , , 

Here, e ^ 0"*" near r = rs{t\ and is zero elsewhere. This has the effect 
of smoothing out the shock into a travelling wave solution. Despite the 
spacetime being pressureless, we can imagine when particles begin to get 
extremely close to one another, as is the case just before a shock begins 
to occur, a force would begin to play a role that kept these particles from 
getting too close. 

Utilizing the viscosity term in the differential equation, we can derive 
a condition on the shock known as the Rankine-Hugionot condition, which 
gives the velocity of the shock (see for e.g. [19]). This is given by 

^ ^ inmt 

dt [M]t 



where 



lim^ — lim . (62) 



This result can be transformed back into the radial coordinate utilized in 
the metric, and the velocity of the shock is found to be 

V,...±^.=l^^-^. (03) 

dt [M]+ ^ ' 

We note that the smoothing either side of the shock necessarily implies that 
the mass is still a monotonically increasing function in the positive radial 
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r 



Figure 1. Initial density profile 



direction. Therefore, rnust necessarily be negative, implying the shock 
travels in the direction towards r = 0. Furthermore, once formed there 
is no mechanism to stop the evolution of the shock, and therefore during 
the evolution it will reach r = 0, at which point the shock will cease to 
exist and the classical solution will be regained. Thus, we see that even if a 
shock forms as a global naked singularity, it will evolve to the end state of 
gravitational collapse as a black hole. 

6. Shock Example 

We wish to highlight the formation and evolution of a shock wave by 
way of example. In particular, specifying an initial matter distribution such 
that the initial tidal force, A(0, r) < for a finite range of r, ensures a shock 
will develop throughout the evolution. We define an initial energy-density 
distribution by the following piecewise continuous function 

f |cos|^-icos^ + l r<r9„ 
[ r>rQg 

While the density profile given (Fig. [T| is an unrealistic starting point 
for gravitational collapse, we note that the energy density in realistic col- 
lapse scenarios may not necessarily be monotonically decreasing, and may 
therefore have corresponding regions where the tidal force takes on the op- 
posite sign. For example, an accretion disc around a black hole will have a 
distribution not dissimilar to the present example. 

While an analytic solution utilizing the above initial conditions can be 
found, it is extremely long, implicit and highly non-linear. However, the 



(64) 
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r 



Figure 2. Characteristic curves of constant mass 




Figure 3. Initial Mass, Mass profile at shock formation, Mass profile after shock forma- 
tion 

classical solution can be plotted even past the initial formation of the shock. 
We do this only to display the nature of the shell crossing singularity being a 
multi- valued function in the mass. Figure [2] displays the characteristic curves 
of constant mass emanating from the initial distribution. The dashed curves 
are those coming from the exterior, vacuum region of the spacetime. The 
shock wave is represented here by a crossing of the characteristic curves. 
Figure |3] displays three plots of the Mass function against the radius of the 
classical solution at three different times. 

7. Conclusion 

A review of the literature associated with LTB collapse shows that major 
phenomenological discussion is focused on SCS (e.g., |H US EOl E] ) . How- 
ever, much of this work is done in terms of arbitrary functions which are 
difficult to interpret. 

The work herein derived LTB regions and Schwarzschild regions in GPG 
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coordinates via an initial value formulation. The advantage of this approach 
is four-fold: 

1. It enabled the metric to be written in terms of a single line element, 
thus avoiding complicated matching schemes at the interface. This dif- 
fers from the standard approach in the literature (e.g., [14') where two 
spacetimes are matched across a boundary by a coordinate transfor- 
mation. Difficulties often arise in showing that this coordinate trans- 
formation is valid everywhere along the boundary as the Jacobian of 
transformation is non-trivial to analyse. By writing the line element 
in a single coordinate patch, no transformation is required and this 
difficulty is avoided. It has already been shown that this method is 
extendable beyond the dust case to include both perfect fluids [9] and 
also a completely general fluid [TO]. However, the equations derived in 
those cases are extremely difficult to handle due to more terms in the 
fluid implying more terms are also required to describe the geometry 
(see [TO] for a discussion of this point). The dust cases analysed in 
this article provides relatively straightforward analytic solutions, and 
hence the structure of the equations is better understood. 

2. Furthermore, as analytic solutions are found from the initial condi- 
tions, this approach is excellent as a test-bed for numerical schemes. 
In particular, the evolution of the boundary of the two regions of the 
spacetime and of the SCS singularity can be tested as exact analytic 
forms of both of these are known. 

3. The metric coefficients are all expressed in terms of the matter fields. 
Furthermore, the dynamics of the spacetime is completely determined 
by two differential equations governing the mass and energy functions. 
The problem is then solved using physically reasonable initial and 
boundary conditions on the energy density. The evolution of this mass 
function can result in multi-valued solutions for particular choices of 
initial conditions, which is exactly equivalent to SCS. The advantage 
of this scheme is that it explicitly identifies the SCS to be equivalent to 
shock waves. This enables one to analyse the dynamical extensions of 
SCS beyond their initial point of formation. In particular, we showed 
that a SCS that forms as a globally naked singularity must become a 
locally naked singularity, and eventually fall to a Schwarzschild singu- 
larity at r = 0, at which point it will cease to effect the system. 

4. The spacetime is easier to visualize due to the physical intuition being 
similar to familiar fluid problems. Rather than dealing with abstract 
four-dimensional coordinates, we deal with the propagation of hyper- 
surfaces, exactly as we imagine the spacetime in which we live. Math- 
ematically, this implies all functions appearing in the solution are not 
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necessarily arbitrary, but are more appropriately based on physical 
intuition. 
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Appendix 



For clarity, we express the extrinsic curvatures Kij, and three-Riemann 
curvatures ^Rij, in terms of the metric coefficients. The Lie derivative oper- 
ator with respect to the normal Cn, is expressed when operating on a scalar 



in equation (13), and we have already let a = 1 (see section 2.1). 



fr -h^nU 










^ 



(Al) 



K 



J.2 ) U^""^^ 



(A2) 



d ff3\ 1 



dr \ r I hi 




(A3) 



3 7? J 



•,dU 
' dr 
















dr 



(A4) 



'2f + '-{U^-U) 
or r ^ ' 



(A5) 



References 



20 



[1] R.J. Adler, J. D. Bjorkon, P. Chen and J.S. Liu, "Simple analytic models of gravi- 
tational collapse", Am. J. Phys. 73 (12) (2005) 1148-59, arXiv:gr-qc/0502040. 

[2] R. Arnowitt, S. Dcscr and C. W. Misncr, "The dynamics of general relativity", in 
Gravitation: An introduction to Current Research (ed. Witten), (John Wiley and 
Sons, Inc., New York, 1962). 

[3] H. Bondi, "Spherically symmetrical models in general relativity", Mon. Not. Roy. 
Astron. Soc. 107 (1947) 410-25. 

[4] D. Christodoulou, "Violation of cosmic censorship in the gravitational collapse of a 
dust cloud", Commun. Math. Phys. 93 (1984) 171-95. 

[5] G. F. R. Ellis, "Relativistic cosmology", in General Relativity and Cosmology (ed. 
B. K. Sachs), (Academic, New York, 1971) 104-82. 

[6] R. Gautreau and J. M. Cohen, "Gravitational collapse in a single coordinate system" , 
Am. J. Phys. 63 (1995) 991-9. 

[7] A. Gullstrand, "AUegemeine losung des statischen einkorper-problems in der ein- 
steinchen gravitations theorie", Arkiv. Mat. Astron. Fys 16 (1922) 1-15. 

[8] H. Muller Zum Hagen, P. Yodzis and H. J. Seifert, "On the occurrence of naked 
singularities in general relativity ii", Commun. Math. Phys. 37 (1974) 29-40. 

[9] P. D. Lasky and A. W. C. Lun, "Generalized lemaitre-tolman-bondi solutions with 
pressure", Phys. Rev. D 74 (2006) 084013. 

[10] P. D. Lasky and A. W. C. Lun, "Spherically symmetric collapse of general fluids", 
Phys. Rev. D 75 (2007) 024031. 

[11] P. D. Lax and B. Wendroff, "Systems of conservation laws", Commun. Pure. Appl. 
Math. 13 (1960) 217-37. 

[12] G. Lemaitre, "L'univers en expansion", Ann. Soc. Sci. Bruxelles A 53 (1933) 51. 

[13] C. W. Misner and D. H. Sharp, "Relativistic equations for adiabatic, spherically 
symmetric gravitational collapse", Phys. Rev. 136 (2) (1964) B571-6. 

[14] C. W. Misner, K. S. Thorne and J. A. Wheeler, Gravitation (Freeman, New York, 
1973). 

[15] R. P. A. C. Newman, "Strengths of naked singularities in tolman-bondi-spacetimes" , 

Class. Quantum Grav. 3 (1986) 527-539. 

[16] B. C. Nolan and F. C. Mena, "Geometry and topology of singularities in spherical 
dust collapse". Class. Quantum Grav. 19 (2002) 2587-605. 

[17] J. R. Oppenheimer and H. Snyder, "On continued gravitational contraction", Phys. 
Rev. 56 (1939) 455-9. 

[18] P. Painleve, "La mccanique classique et la theorie de la relativite" , C. R. Acad. Sci. 
173 (1921) 677-80. 

[19] J. SmoUer, Shock Waves and Reaction-Diffusion Equations (Springer- Verlag New 
York Inc., New York, 1983). 



21 



[20] P. Szckcros and A. Lun, "What is a shell crossing singularity", J. Austral. Math. 
Soc. Ser. B. 41 (1999) 167-79. 

[21] R. C. Tolman, "Effect of inhomogeneity on cosmological models", Proc. Nat. Acad. 
Set. USA 20 (1934) 169-76. 

[22] R. Wald, General Relativity (U. of Chicago Press, 1984). 

[23] P. Yodzis, If. J. Seifert and H. Muller zum Hagen, "On the occurrence of naked 
singularities in general relativity", Commun. Math. Phys. 34 (1973) 135-48. 



22 



